Gaussian Processes (GP) have established themselves as a powerful and flexible tool for modelling nonlinear functions within a Bayesian non-parametric framework [Williams and Rasmussen, 2006]. Multi-output Gaussian processes (MOGP) generalise this powerful framework to the vector-valued random field [Alvarez et al., 2012] by capturing correlations not only across different inputs but also across different output functions. This characteristic has been experimentally shown to provide better predictions in fields such as geostatistics [Wackernagel, 2003], heterogeneous regression [Moreno-Muñoz et al., 2018], and the modelling of aggregated [Yousefi et al., 2019] and hierarchical datasets [Ma et al., 2023]. The primary focus in the literature on MOGP has been on developing an appropriate cross-covariance function between the multiple outputs. Two classical approaches for defining such cross-covariance functions are the Linear Model of Coregionalization (LMC) [Journel and Huijbregts, 1976] and process convolutions [Higdon, 2002]. In the former, each output corresponds to a weighted sum of shared latent random functions.